rm(list = ls())
getwd()
setwd("C:/Users/griff/OneDrive/Desktop")

library(tidyr)
library(dplyr)
library(ggplot2)

#######FOR FIGURE 1.1#######

#some intial data wrangling with the full GTD and V-Dem data in time series format
data <- read.csv(file="Final Data and R Scripts/ISA_ts_data.csv", na.string = "")
vdem <- read.csv(file="Final Data and R Scripts/V-Dem v10/othervdem.csv", na.string = "")
merge1 <- merge(vdem, data, by=c("Country","Year"), all.x = TRUE)
merge1 <- merge1 %>% na_if("NA")

merge1$Insurgency <- as.numeric(merge1$intensity >=1)
merge1$Insrgency[merge1$Insurgency==1] <- "Insurgency"
merge1$Insrgency[merge1$Insurgency==0] <- "No Insurgency"
aggregate (merge1$SQ_attacks, by=list(Category=merge1$Insurgency), FUN=sum)

merge1 <- merge1 %>% dplyr::na_if(-66)
merge1 <- merge1 %>% dplyr::na_if(-77)
merge1 <- merge1 %>% dplyr::na_if(-88)

#Creating Figure 1.1
ggplot(data=subset(merge1, !is.na(Insrgency))) + 
  geom_bar(aes(x=Insrgency, y=SQ_attacks), stat = "identity") +
  labs(title = "Figure 1.1", subtitle = "Status-Quo Terrorism and Insurgency", 
       x = "", y = "Status-Quo Attacks") +
  theme_minimal()+
  theme(plot.title = element_text(hjust = 0.5), plot.subtitle = element_text(hjust = 0.5))


#######FOR FIGURE 1.2########
data1 <- read.csv(file="Liberal Constraints/JTPV/Data/SQ_attacks_during_insurgency.csv", na.string = "")
str(data1)
data1$IGO_Context <- as.numeric(as.character(data1$IGO_Context))
data1$loggdp <- as.numeric(as.character(data1$loggdp))
data1$polity <- as.numeric(as.character(data1$polity))
data1$Polity <- as.numeric(as.character(data1$Polity))

require(pscl)
require(MASS)
require(boot)

n4 <- zeroinfl(SQ_attacks ~ Polity + IGO_Context + softpta +
                 territorial_control + loggdp + past_incidents,
               data=data1, dist="negbin")


pred <- data.frame(Polity = rep(seq(from = min(data1$Polity, na.rm=TRUE), to = max(data1$Polity, na.rm=TRUE), length.out = 100)),
                   IGO_Context = mean(data1$IGO_Context, na.rm=TRUE), 
                   softpta = mean(data1$softpta, na.rm=TRUE),
                   territorial_control = mean(data1$territorial_control, na.rm=TRUE),
                   loggdp = mean(data1$loggdp, na.rm=TRUE),
                   past_incidents = mean(data1$past_incidents, na.rm=TRUE))

pred$SQ.pred <- predict(n4, pred, type = "response")
pred

ggplot(pred, aes(x = Polity, y = SQ.pred)) +
  geom_smooth(mapping = aes(x = Polity, y = SQ.pred), color = "blue", size = 1) +
  labs(title = "Figure 1.2", subtitle = "Predicted Status-Quo Count by Polity Score",
       x = "Polity Score", y = "Predicted Status-Quo Attacks") +
  theme_classic() +
  theme(plot.title = element_text(hjust = 0.5), plot.subtitle = element_text(hjust = 0.5))
